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Chapter  1 
Introduction 


Damage  due  to  an  external  impact  is  a  serious  problem  for  fiber  reinforced  organic 
matrix  composite  structures.  An  impact  may  cause  delamination  or  other  damage  to 
the  structures.  The  structures  may  require  routine  inspections  to  insure  its  integrity. 
A  system  that  could  automatically  detect  and  report  an  impact  would  make  inspection 
for  impact  damage  only  necessary  when  an  impact  of  significant  energy  is  received. 

Several  authors  have  published  methods  to  identify  an  impact  on  a  plate  using 
sensor  measurements.  Wu  et.  al.  [1],  [2]  and  [3],  presented  an  optimization  procedure 
to  calculate  the  impact  history  of  an  impact  at  a  known  location  on  a  composite 
plate.  They  used  an  array  of  strain  gauges  to  measure  the  plate  response.  They 
also  developed  a  method  to  select  the  impact  location  from  several  possibilities  on  an 
aluminum  plate.  However,  they  did  not  develop  a  complete  system  for  locating  and 
reconstructing  an  unknown  impact  on  a  composite  plate. 

Another  way  of  solving  the  identification  problem  is  using  a  neural  network. 
Shaw  [4]  and  Jones  [5]  developed  a  neural  network  based  identification  method  and 
successfully  found  the  location  of  a  known  force.  The  training  of  the  neural  network 
becomes  impractical  when  the  force  history  is  unknown.  Reconstructing  an  impact 
force  history  would  require  nearly  an  infinite  number  of  training  impacts. 

Approaching  the  problem  of  impact  identification  in  another  way,  Gunther  et 
al  [6],  recorded  the  time  of  arrival  of  the  signal  at  four  different  fiber-optic  sensors 
and  triangulated  to  find  the  non-damaging  impact  location.  However,  in  practice. 
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using  details  of  the  signal’s  time  of  arrival  is  difficult.  Because  the  wave  is  dispersive, 
the  initial  part  of  the  sensor  response  has  a  low  amplitude  and  a  high  frequency 
content.  If  there  is  noise  in  the  sensor,  it  is  difficult  to  differentiate  between  the 
arrival  of  the  signal  and  the  sensor  noise. 

A  robust  and  accurate  impact  identification  system  for  a  beam  was  developed  by 
Choi  and  Chang  [7].  An  optimization  technique,  a  smoother /filter,  was  used  to  solve 
the  inverse  problem. 

The  smoother /filter  is  employed  in  this  research.  Its  use  is  expanded  to  the  two- 
dimensional  impact  identification  problem  on  a  composite  plate.  This  approach  will 
provide  a  robust  and  reUable  means  of  identifying  the  impact  on  composite  plates 
with  clamped  or  free  boundary  conditions  and  noisy  sensor  data. 


Chapter  2 

Problem  Statement 


The  problem  is  identification  of  impacts  on  composite  plates  using  a  distributed 
sensor  array.  The  identification  includes  determination  of  the  impact  location  and 
reconstruction  of  its  force-time  history.  The  duration  of  the  low-mass,  high-speed 
impact  is  on  the  order  of  one  millisecond.  The  composite  plates  are  general  symmetric 
laminates.  The  sensors  measure  the  plate’s  response  to  an  impact. 

Specifically,  the  research  addresses  the  identification  problem  for  the  following 
conditions: 

•  The  composite  plate  has  a  symmetric  layup. 

•  The  composite  plate  has  clamped  or  free  boundary  conditions  but  does  not  have 
any  stiffeners. 

•  The  length  and  width  of  the  composite  plate  are  much  greater  than  its  thickness 
so  that  shear  effects  on  plate  bending  are  negligible. 

.  •  The  impact  does  not  cause  damage  to  the  composite. 

•  The  impact  occurs  at  a  point.  The  size  of  the  impactor  is  much  less  than  the 
dimensions  of  the  plate. 

•  The  sensors  are  small  and  widely  spaced  so  that  the  bending  stiffness  of  the 
plate  is  not  effected. 
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The  identification  system  is  introduced  in  Figure  1.  Using  the  measured  response 
of  the  plate  and  the  application  parameters,  the  identification  system  reports  the 
impact  location  and  the  reconstructed  force  time  history. 
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Figure  1:  The  impact  identification  system.  The  system  identifies  the  un¬ 
known  impact  on  the  composite  plate.  With  the  plate  response 
and  application  parameters,  the  system  reports  the  impact  loca¬ 
tion  and  force  time  history. 


Chapter  3 
Approach 


As  described  in  the  problem  statement,  the  identification  system  identifies  the 
unknown  impact  on  a  composite  plate.  The  input  to  the  system  is  the  plate  response 
and  the  application  parameters  which  describe  the  plate  and  sensor  properties.  The 
output  is  the  impact  identification  which  includes  the  impact  location  and  force  time 
history. 

This  research  includes  the  theoretical  development,  the  computer  implementation 
and  the  experimental  testing  of  the  identification  system. 

The  theory  was  developed  for  the  four  components  of  the  identification  system,  see 
Figure  2.  The  sensors,  the  preprocessor,  the  forward  model  and  the  inverse  problem 
solver. 

•  The  plate  response  is  measured  using  small  circular  piezoelectric  sensors.  The 
sensors  measure  strain  where  bonded  to  the  surface  of  the  composite  plate. 

•  The  preprocessor  provides  experimental  data  to  the  inverse  problem  solver.  The 
preprocessor  determines  which  sensors  are  closest  to  the  impact.  These  sensors 
are  isolated  so  that  the  most  relevant  data  is  used  in  the  solution  process. 

•  The  forward  model  calculates  the  plate  response  as  a  bending  wave  radiating 
from  the  point  of  impact.  The  bending  wave  is  modeled  using  Kirchoff ’s  theory 
for  composite  plates.  The  governing  equation  is  solved  using  a  discrete  time 
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integration.  The  model  response  is  compared  to  the  experimental  measured 
response  in  the  inverse  problem  solver. 

•  The  inverse  problem  solver  identifies  an  impact  that  will  minimize  the  difference 
between  the  model  response  and  the  measured  response.  The  minimization  is 
accomplished  with  an  outer  loop  estimating  the  impact  location  and  an  inner 
loop  reconstructing  the  force  history,  see  Figure  3.  The  inner  loop  provides  a 
figure-of-merit  that  is  used  to  update  the  estimated  location.  The  procedure 
continues  until  the  the  best  estimate  of  the  impact  location  is  found. 

The  identification  system  was  implemented  with  the  computer  code  IDIMPACT. 
The  system  runs  on  one  machine  that  collects  and  digitizes  the  data  and  then  cal¬ 
culates  the  solution.  The  code  solves  the  identification  problem  given  the  digitized 
response  of  the  composite  plate.  It  reports,  in  a  graphic  interface,  the  impact  location 
and  force  time  history. 

The  ability  of  the  identification  program,  IDIMPACT,  to  identify  an  unknown 
force  was  demonstrated  on  a  composite  plate.  The  system  was  then  extensively 
tested  with  a  distributed  array  of  impacts  to  determine  its  accuracy  and  reliability. 
Finally,  the  robustness  of  the  identification  system  was  tested  by  introducing  noise  in 
the  sensor  data. 
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Identification  system 


Figure  2:  Inside  the  identification  system.  The  sensors  measure  the  plate 
response,  the  preprocessor  isolates  the  relevant  sensor  region,  and 
the  inverse  problem  solver  identifies  the  impact  using  the  measured 
response  and  the  model  response. 
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Figure  3:  The  inner  and  outer  loop  of  the  inverse  problem  solution. 


Chapter  4 

Sensor  and  Preprocessor 

The  sensors  and  the  preprocessor  provide  experimental  data  to  the  inverse  problem 
solver.  Piezoelectric  sensors  were  selected  to  measure  the  response  of  a  composite 
plate  to  an  impact.  The  response  is.digitized  by  the  analog  to  digital  converters.  The 
preprocessor  scans  the  data  to  find  the  sensors  closest  to  the  impact.  The  closest 
sensors  and  the  corresponding  area  are  used  in  the  solution  process. 


4.1  Piezoelectric  Sensors 


The  small  circular  piezoelectrics  are  manufactured  by  Piezo  Kinetics  Incorporated. 
The  sensors  are  good  for  dynamic  reading  because  they  have  a  very  fast  response  time 
and  have  sufficient  voltage  output  that  sophisticated  amplification  is  not  necessary. 

The  relation  between  the  voltage  output  and  the  strain  is  supplied  in  the  manu¬ 
facturer’s  application  notes.  The  one  dimensional  strain  relation  is 


€xx  — 


Al 

I 


(4.1) 


The  strain  Cxx  is  the  ratio  of  i,  the  length  of  the  piezo,  over  Al,  the  change  in 
length.  V  is  the  voltage  and  tpz  is  the  thickness  of  the  piezo.  The  piezoelectric 
constant,  dai,  relates  the  mechanical  strain  (fractional  deformation  of  the  ceramic)  to 
the  electrical  voltage.  This  constant  is  supplied  by  the  manufacturer.  The  piezo  is 
poled  perpendicular  to  the  x-y  plane.  See  Figure  4  for  an  explanation  of  the  geometry. 
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In  this  research,  small  circular  piezoelectrics  are  used.  The  diameter,  d,  of  the 
piezo  is  .25  inches.  The  thickness,  t,  is  .01  inches.  The  equation  to  obtain  voltage 
from  a  two  dimensional  strain  field  is 

V  =  JT  i^xx  +  Cyy)  (4-2) 

“31  (1  — 

This  equation  simplifies  to  the  one-dimensional  equation  when  the  stress  ay  is  set 
to  zero.  In  that  case,  €yy  =  — i^Cn,  and  Eq.  4.2  becomes 

V  =  (4.3) 

“31 

The  experimental  strain  is  attained  from  the  experimental  voltage  data  using 

Seip  —  Se  ^  ^ exp  (4-4) 

where  Vexp  is  the  voltage  data  from  each  sensor,  is  the  gain  to  convert  the  voltage 

to  strain  and  includes  the  constants  of  Eq.  4.2  and  the  effects  of  bonding.  This 

constant  can  be  determined  experimentally  by  comparing  the  forward  model  to  the 
experimental  data.  The  experimental  strain,  Seap,  contains  the  experimental  strain 
data  for  all  the  sensors. 

The  circular  piezos  measure  a  non-directional  strain  property.  The  combination 
of  strain,  Cxx  +  (-yy,  is  a  two-dimensional  strain  invariant.  The  measurement  is  inde¬ 
pendent  of  the  orientation  of  the  coordinate  system.  This  can  be  demonstrated  using 
the  two-dimensional  strain  transformation  tensor  where  6  is  the  angle  of  coordinate 
rotation. 

cos^  (^)  sin^  (0)  cos  {6)  sin  (^)  txx 

sin^  {6)  cos^  {9)  -  cos  {9)  sin  {9)  Cyy  (4.5) 

—2  cos  {9)  sin  (9)  2  cos  (0)  sin  (9)  cos^  (9)  —  sin^  (9)  _  ^xy  _ 

The  strain  combination  is  then 


^xx  ^yy  ^xx  4“  Cyy 


(4.6) 
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The  quantity  is  not  changing  with  the  rotation  of  the  coordinate  system.  By  providing 
a  non-directional  means  of  measuring  the  plate  response,  the  sensor  simpUfies  the 
inverse  problem  solution. 

The  sensors  used  in  this  research  measure  the  strain  invariant  +  Cyy.  The 
experimentally  measured  quantity  becomes 

Sgajp  ”  ^xx  ^yy 

It  is  assumed  that  the  diameter  of  the  circular  piezo  is  significantly  smaller  than 
the  wave  length  of  the  plate  response.  The  sensor  response  is  calculated  as  a  point 
property  of  strain.  It  is  assumed  the  efiect  of  averaging  the  experimental  strain  over 
the  sensor  surface  area  is  small.  This  assumption  will  be  tested  in  the  application  of 
the  identification  system  by  comparing  the  sensor  diameter  to  the  global  response  of 
the  composite  plate. 

4.2  Preprocessor 

The  preprocessor  selects  the  sensors  closest  to  the  impact  location.  This  allows 

/ 

for  a  region  of  the  sensor  array  to  be  isolated  for  the  inverse  problem  solver. 

The  preprocessor  scans  the  experimental  strain  to  find  the  sensors  closest  to  the 
impact.  It  determines  the  close  sensors  by  comparing  the  energy  in  an  initial  time 
window  of  the  data.  It  calculates  the  energy  using 

E  =  e'sLpCO  0-8) 

t=0 

where  E  is  the  energy  of  the  sensors  being  scanned.  Ng  is  the  number  of  data  points 
included  in  the  initial  time  window.  The  number  of  points  used  should  be  enough 
that  sensors  close  to  the  impact  have  a  large  energy  content  while  those  far  away  have 
close  to  zero, 

The  energy  of  each  of  the  sensors  are  compared.  The  sensors  with  the  greatest 
energy  and  .the  corresponding  region  of  the  plate  are  then  used  by  the  inverse  prob¬ 
lem  solver.  This  method  provides  a  method  that  is  not  susceptible  to  noise  in  the 
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experimental  data  and  provides  a  reliable  way  to  isolate  a  region  of  the  sensor  array. 
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Figure  4:  The  description  of  the  relationship  between  strain  and  voltage  for 
the  piezoelectric.  The  piezo  is  poled,  P,  in  the  vertical  direction 
and  strained  horizontally. 


Chapter  5 

Forward  Model 


The  forward  model  calculates  the  plate  response  as  a  bending  wave  radiating 
from  the  point  of  impact.  The  bending  wave  is  modeled  using  Kirchoff’s  theory 
for  composite  plates.  A  governing,  equation  is  developed  that  describes  the  plate 
response,  both  displacement  and  strain,  at  specific  points.  The  governing  fourth- 
order  differential  equation  is  reduced  to  a  system  of  first-order  equations  using  a 
Fourier  transformation.  The  equations  are  recast  in  state-space  form  and  solved  with 
a  discrete-time  integration.  The  plate  can  have  either  clamped  or  free  boundary 
conditions.  The  model  response  is  compared  to  the  experimental  measured  response 
in  the  inverse  problem  solver. 

The  accuracy  of  the  forward  model  is  verified  by  comparison  with  an  analytic 
solution  of  an  isotropic  plate  and  with  an  experimental  impact  on  a  composite  plate. 
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5.1  Governing  Equation 

The  governing  equation  for  the  composite  plate  can  be  expressed  in  state  space 
form.  Starting  with  Kirchoff’s  plate-bending  theory  generalized  for  composite  plates 
with  symmetric  layups,  the  governing  equation  for  a  point  impact  is 


+  2(Di2  +  2D66) 


,  d^w  _  d^w 

^  d^w  _  d^w  _  d^w 

+  4T)26  t  +  D22-:r-r  +  771 


d*w 


(5.1) 


’  dxdy^ 


dy^ 


dt^ 


f(t)S{x,y)  is  the  point  load  located  at  (0, 0),  Dij  are  elements  of  the  composite  lam¬ 
inate  matrix  D,  (Tsai  [8]),  and  w  is  the  out-of-plate  deflection  of  the  plate. 

This  governing  equation  is  applied  to  an  infinite  plate  by  using  the  two-dimensional 
discrete  Fourier  transformation  with  respect  to  the  spatial  variables  x  and  y. 


w{xi,yj,t) 


iV-1  M-1 

n=0  m=0 


(5.2) 


The  frequency  components  are  calculated  using  Km  —  and  /Cn  =  ^-  X  and 
Y  'are  the  extent  of  the  spatial  domain.  Wmn  is  the  amplitude  of  each  frequency  com¬ 
ponent.  Applying  the  transform  to  both  sides  of  the  equation  yields  the  component 
equation 


{DnK%^  -t-  ADi^K^Kn  +  2(I?i2  +  2D^)k^k\+ 
4D26l^ml^n  +  ^22l^n)  ~ 


(5.3) 


These  component  equations  can  be  regrouped  by  introducing  the  state-space  formu¬ 
lation. 


d  ‘lOmn 
^  ihmn 


0  1 
—Kmn  0 


'^mn 

Wmn 


(5.4) 


where 
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uj{n  +  Nm)  =  uJmn  (5-9) 

and  K  is  a  diagonal  matrix,  [N  •  M  x  N  •  M],  with  the  diagonal  elements 

K(n  +  Nm,  n  +  Nm)  =  Kmn  (5.10) 

where  m  =  0,  . . . ,  (M  —  1)  and  n  =  0,  . . . ,  (A^  —  1) 


Chapter  5.  Forward  Model 


18 


5.2  Plate  Response 

A  state  space  formulation  is  used  to  calculate  the  plate  response.  The  Fourier 
Transform  provides  the  out-of-plane  deflection  of  the  composite  plate  at  a:  =  Xj  and 

y  =  Vj 

N-lM-l 

w{Xi,yj,t)=  (5-11) 

n=0  m=0 

The  equation  can  be  represented  in  matrix  form.  u>  is  a  [1  x  N  •  M]  vector.  A 
vector  with  [N  •  M  xl]  elements  is  defined  for  the  exponential  terms.  For  m  = 
0,  . . . ,  (M  -  1)  and  n  =  0,  . . . ,  '{N  -1) 

C^(rH- Nm)  =  (5.12) 

The  calculation  of  the  displacement  becomes 

w{xi,yi,t)  =  Cy,Q  (5.13) 

A  more  general  form  of  this  equation  can  be  defined 

f 

s  =  Cz  (5.14) 

where  s  is  the  general  measured  quantity  -  displacement,  acceleration  or  strain,  s  is 
a  [NS  X  1]  vector  where  NS  is  the  number  of  sensors  in  the  system.  Each  element 
of  s  is  the  response  of  one  sensor  at  time  step  n.  The  C  matrix,  is  the  observation 
matrix.  For  one  sensor 


C=[c„  o]  (5.15) 

C  is  a  [NS  X  N  ■  M]  matrix  where  each  row  represents  one  sensor. 

The  observation  matrix  can  be  simplified  using  the  symmetry  of  the  impact  prob¬ 
lem  and  the  properties  of  a  Fourier  transformation  for  real  and  even  functions,  see 
Appendix  A.  The  simplified  submatrix  becomes  a  ^NS  x  (y  1)  •  (y  +  1)]  matrix. 
For  m  =  0,  . . . ,  y  and  n  =  0,  . . . ,  y  ,  each  row  becomes 
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[  1  •••  2cos(/c„j/j) 

2  cos(KmXi)  •  •  •  4  cos(«m3^i)  cos(/c„2/j) 

COS(KM/2^i)  2cOS(/CM/2a:i)cOS(K„yj) 


cos(Kiv/2yj) 

2  COs(KmXi)  COs{KN/2yj) 

C0s{KM/2Xi)  COs(/Cjv/2yi)  ] 

(5.16) 


5.3  strain  Measurement 


Strain  is  often  the  measured  quantity  in  structures  applications.  To  calculate 
strain  with  the  measured  quantity  s{t),  an  appropriate  strain  observation  matrix  is 
required. 

The  equations  for  surface  strain  are 

h 


2  dx^ 

^  n=0  m-0 


—  11  V'  1^2  ^7i  pi«mXj+i/Cn2/j 

“  o  2^  ^m^rnn^ 


w 


2  ^2/2 

=  9  E  E  «nWmn 

^  n=0  m=:0 


gt/CmXi+i/Cnyj 


(5.17) 


where  h  is  the  thickness  of  the  plate,  Cu  is  the  strain  measured  on  the  surface  of  the 
plate  in  the  x  direction,  and  tyy  is  the  strain  in  the  y  direction. 

Similar  to  displacement,  the  strain  equations  can  be  written  in  matrix  form.  The 
strain  equations  become 


^xx(.Xi,yj,t')  —  Cxx^ 


(5.18) 
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^yy^  (5.19) 

The  elements  of  the  the  strain  submatrices  Cxx  and  Cyy  are 

{Cxx){n  +  Nm)  =  K^(CJ(n  +  iV7n)  , 

(5.20) 

{Cyy){n  +  Nm)  =  Kl{Cy;){n  + Nm) 

where  m  =  0,  . . . ,  {M  —  1)  and  n  =  0,  . . . ,  {N  —  1)  ■ 

The  general  measured  quantity  s(t)  can  now  be  used  to  calculate  the  strain.  Al¬ 
though  the  observation  equation  remains  unchanged, 

s(t)  =  Cz  (5.21) 

the  observation  submatrix  has  changed.  For  each  sensor  measuring  strain  in  the  x 
direction  and 

C=[c.x  o]  (5.22) 

and  for  each  sensor  measuring  in  the  y  direction 

■'  C=[c,,  o]  (5.23) 

5.4  Discrete  Time  Domain 

To  solve  the  governing  equation,  the  state-space  formulation  developed  in  the 
previous  section  is  converted  to  the  discrete  time  domain.  The  system  formulation  in 
the  continuous  time  domain  is 


z  =  Az  -f-  B/(t)  (5-24) 


s(t)  =  Cz 


(5.25) 
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z(n  +  l)  = 


governing  equations  become 

<f>ll  <f>12 

z(n)+ 

7i 

021  022 

72 

f{n) 


(5.26) 


=  $z(n)  +  r/(n) 

and  the  observation  equation  becomes 

s(n)  =  Cz(n)  (5.27) 

where  n=  0,  . . . ,  {NT  —  1)  •  The  new  matrices  are  defined  as 

^  =  exp  (ATj)  (5.28) 

r  =  /  exp(At)dtB  (5.29) 

Jo 

/ 

See  Franklin,  et.  al.  [9]  for  a  discussion  of  these  relations.  $  is  a  block  diagonal 
matrix,  through  <^22  are  [N  •  M  x  N  •  M]  diagonal  matrices.  71  and  71  are 
[iV  •  M  X  1]  vectors.  The  observation  matrix  C  remains  unchanged. 

The  solution  to  the  discrete  time  equations  begins  with  an  initial  condition,  z(0) 
and  the  initial  value  of  the  force  /(O).  The  variables  are  calculated  at  the  discrete 
time  step  n  where  the  time  is  t  =  nTg.  Tg  is  the  sampling  period.  If  the  plate  is 
initially  at  rest,  z(0)  =  0.  The  time  history  of  the  strain  is  constructed  by  calculating 
z{n  +  1)  from  z(n)  and  /(n)  and  updating  the  measured  quantity  s(n). 

The  discrete  time  relations,  Eq.  5.28  and  Eq.  5.29,  can  be  solved  to  facilitate  the 
numeric  implementation.  Defining  A^ut,  a  2  x  2  subarray  of  the  A  matrix, 

0  1 
— K(n  +  Nm,  n  +  Nm)  0 


Ajub  — 


(5.30) 
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a  subarray  ^sub  of  the  $  matrix  can  be  calculated  using  Eq.  5.28.  For  n  and  m^O, 
the  resulting  $  subarray  is 


^sub 


<f>ii{n  +  Nm)  4>i2in +  Nm) 
(l>2i{n  +  Nm)  <f)22{n  +  Nm) 

cos  (Tgy/Kjfir^ 
y/KmnSin  (TsVKmn) 


1 

y/^mn 

COS 


sift  (Ts^/ Kfnn^ 

(t.vk;;;;) 


(5.31) 


Similarly,  defining  Bsub  as  a  subarray  of  the  B  matrix,  the  elements  of  the  F,  Tsut, 
can  be  calculated  using  Eq.  5.29.  The  B  subarray  is 


and,  for  n  and  m  /  0,  the  resulting  F  subarray  is 


(5.32) 


F  sub 


4>i{n  +  Nm) 
(t>2{n  +  Nm) 


cos  {TssJ Kfjir^  d" 

■l^mn  -^rt 

.  sin  (Ts'^ Kjjir^ 

V  ^mn 


If  m  and  n  =  0,  then 


^sttb  — 


F  sub  — 


1  Ts 
0  1 


bTs 


(5.33) 


(5.34) 


(5.35) 
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5.5  Forward  Model  Verification 

To  verify  the  accuracy  of  the  forward  model  response,  it  was  compared  to  an 
analytic  solution  of  an  isotropic  plate  and  the  measured  response  of  a  composite 
plate. 

The  analytical  solution  for  an  isotropic  aluminum  plate  is  presented  in  Doyle  [10]. 
The  strain  response:  is  calculated  in  radial  coordinates  using  the  Fourier  Transform 

=  (5.36) 

j=0 

where 

i„(j)  =  [(•'o  -h)-i  (n  -  1^2)  +  f  -  -^2)]  (5.37) 

and 

Fij)=j:f{p)e=^  (5.38) 

p=0 

/  is  the  force  time  history.  D  is  the  bending  stiffness  of  the  isotropic  plate,  i  is  equal 
to  The  argument  for  the  Bessel  functions  J,  Y  and  K,  is 

z  =  Pr  0  =  y/uj  ^  (5.39) 

where  p  is  the  density  of  the  isotropic  material. 

A  half-sine-wave  impact  was  the  input  to  the  forward  model  and  the  analytic 
solution.  The  responses  are  compared  in  Figure  5.  A  very  small  error  was  seen  to 
accrue  in  the  forward  model. 

The  forward  model  was  also  compared  to  the  measured  response  from  an  actual 
impact  on  composite  plate.  A  T800  composite  plate  with  a  [45/90/  -  45/O2/45/O2/  - 
45/0]i  layup  and  clamped  boundaries  was  used.  The  composite  D  matrix  for  the 
plate  is 

2265.4  641.6  185.1 

641.6  1990.2  185.1 

185.1  185.1  760.6 


(5.40) 


Chapter  5.  Forward  Model 


24 


where  the  units  are  Ibs-in.  The  mass  per  unit  area  of  the  plate,  fh  is  .000165  lb  s^/ft^. 
The  plate  is  30”  x  36” .  The  thickness  of  each  ply  is  .007  inches  for  a  total  thickness 
is  .14  inches. 

The  response  was  calculated  with  Eq.  5.26  and  Eq.  5.27  where  N  =■  M  =  180, 
NT  =  45,  Ts  =  4.48  x  10“®  seconds,  and  X  =  T  =  90  inches. 

The  plate  was  hit  at  its  center  and  the  strain  response  was  measured  at  four 
locations.  The  plate  dimensions,  sensors  coordinates  and  boundaries  are  shown  in 
Figure  6.  The  sensors  were  located  along  the  45  degree  ply  direction  and  the  90 
degree  ply  direction.  The  sensors  at  different  ply  directions  test  if  the  model  includes 
the  anisotropic  properties  of  the  plate.  - 

The  impact  force  was  recorded  using  an  instrumented  hammer.  The  hammer  was 
manufactured  by  PCB  Piezotronics.  It  has  a  208  A04  force  transducer  and  a  480A 
power  unit.  The  impact  force  history  is  shown  in  Figure  7.  The  model  response  to 
this  impact  is  compared  to  the  experimental  response  in  Figure  8  and  Figure  9.  The 
forward  model  was  shown  to  be  a  simple  and  reasonable  means  of  calculating  the 
impact  response  of  the  plate. 

/ 

5.6  Boundaries 

In  real-Ufe  applications  the  plate  will  have  boundaries.  The  effects  of  the  bound¬ 
aries  can  be  approximated  by  adding  mirror  images  of  the  actual  impact  about  the 
boundaries. 

A  sensor  bonded  to  an  actual  plate  with  boundaries  measures  the  initial  wave 
generated  by  the  impact  and  the  wave  after  it  has  reflected  from  the  boundaries.  In 
the  forward  model,  the  effect  of  these  boundaries  can  be  simulated  by  adding  the 
infinite  plate  solution  with  the  impact  located  symmetrically  about  the  boundaries. 
The  sensor  in  the  forward  model  then  includes  the  response  from  the  impact  as  well  as 
the  response  from  the  image  impacts.  The  sum  of  the  impact  responses  approximate 
the  effect  of  the  plate  boundary. 

The  location  of  the  image  impacts  to  simulate  the  first  reflection  from  each  of 
the  four  boundaries  is  shown  in  Figure  10.  The  location  of  the  boundary  is  measured 
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Milliseconds 


Figure  5:  Comparison  of  the  forward  model  and  the  analytic  solution.  The 
*■  response  of  an  isotropic  plate  to  a  half  sine  wave  impact  is  cal¬ 
culated  using  the  forward  model  and  the  analytic  solution.  The 
location  of  the  sensor  relative  to  the  impact  is  shown.  A  very 
small  error  accrues  in  the  forward  model. 
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Figure  6;  The  experimental  setup  to  verify  the  forward  model.  The  zero 
direction  is  in  the  x  direction. 
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Figure  7;  Experimentally  recorded  impact.  This  impact  was  used  to  compare 
the  forward  model  response  and  the  experimental  response.  The 
x-axis  is  the  zero  ply  direction. 
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Sensor  at  x=.71,  y=.71 


Milliseconds 


Sensor  at  x=1.4,  y=1.4 


Milliseconds 


Figure  8:  Verification  of  the  forward  model.  The  graphs  compare  the  model 
response  to  the  experimental  response  at  two  locations. 


(9OIX)  (9OIX)  ^ 
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Sensor  at  x=0,  y=1.0 


Sensor  at  x=0,  y=2.0 


Milliseconds 


Figure  9;  Verification  of  Composite  Plate  Model 
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perpendicular  from  the  impact  location.  Image  impacts  can  also  be  added  to  simulate 
second  and  third  reflections. 

In  the  forward  model,  only  the  observation  matrix  is  affected.  The  C  matrix  is 
the  summation  of  the  observation  matrices  translated  to  the  correct  point  relative  to 
each  image  impact.  The  observation  matrix  for  the  first  reflection  from  the  clamped 
boundaries  is 


C  —  C(xs,  2/s)  +  C(a;s  2x(,„,2/s  ^Vbn)  (^•‘^1) 

where  {xb^,ybj  is  the  location  of  the  boundary  measured  perpendicular  from  the  im¬ 
pact  location.  There  are  a  set  of  coordinates  for  each  boundary,  so  n  =  0,  . . . ,  {Nb  -  1) 
and  Nb  is  the  number  of  boundaries.  For  free  boundary  conditions 

C  =  C(xs,  Vs)  -  C(xs  -  2xb„,ys  -  22/6„) 

where  the  only  difference  is  the  sign  of  the  image  impact. 

/The  forward  model  with  the  image  impact  is  compared  to  the  measured  response 
in  Figure  11.  The  impact  and  sensor  were  located  so  there  was  significant  boundary 
effect  in  the  signal.  The  impact  was  near  the  center  of  the  plate,  12.5  inches  from 
right  clamped  boundary  and  20  inches  from  the  top.  The  sensor  was  located  5  inches 
from  the  right  boundary  and  25  inches  form  the  top.  With  the  sensor  further  away 
from  the  impact,  the  match  between  the  model  and  the  measured  response  is  not  as 
close  as  in  Figure  8  and  Figure  9.  The  general  shape  of  the  response  is  present  but 
some  of  the  higher  frequency  reflected  waves  are  not  well  modeled. 


5.7  Global  Response 

The  impact  response  of  the  entire  plate  can  be  calculated  at  one  point  in  time. 
At  every  point  in  the  space  domain  Eq.  5.21  is  used  to  calculate  the  plate  deflection. 
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f(2y53,2yb3)  | 

I 


Figure  10:  Simulation  of  clamped  boundaries.  The  boundaries  of  the  real 
plate  are  simulated  using  image  impacts.  The  image  and  the 
actual  impact  are  located  symmetrically  about  the  boundaries. 
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Milliseconds 


Figure  11:  Verification  of  the  forward  model  where  there  is  significant  bound¬ 
ary  effect.  The  sensor  is  5  inches  from  the  right  clamped  boundary 
and  25  inches  form  the  top.  The  impact  is  12.5  inches  from  the 
right  boundary  and  20  inches  from  the  top. 
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The  force  time  history  of  Figure  7  was  applied  to  the  center  of  the  composite  plate. 
The  response  at  time  step  n  =  22  or  t  =  .986  milliseconds  is  shown  in  Figure  12. 

From  the  global  response,  it  is  evident  that  the  wavelength  of  the  response  is 
several  inches.  With  this  wave  length,  the  assumption  of  the  sensors  measuring  a 
point  property  of  the  strain  appears  reasonable. 
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Figure  12:  The  global  response  of  the  T800  composite  plate.  The  plate  with 
clamped  boundaries  is  hit  in  the  middle.  The  displacement  is 
calculated  at  time  step  n  =  22  or  t  =  .986  milliseconds. 


Chapter  6 

Inverse  Problem  Solver 


The  inverse  problem  solver  is  the  key  component  of  the  identification  system. 
The  solver  identifies  an  impact  that  will  minimize  the  difference  between  the  model 
response  and  the  measured  response. 

The  solution  to  the  inverse  problem  has  two  major  elements.  The  outer  loop 
searches  for  the  impact  location  and  the  inner  loop,  the  smoother/filter,  reconstructs 
the  impact  force  time  history.  The  solution  starts  from  the  outer  loop  with  an  estimate 
of  the  impact  location.  At  this  location  the  filter /smoother  reconstructs  a  force 
history.  The  solution  continues  when  the  figure-of- merit,  based  on  the  difference 
between  the  model  response  and  the  experimental  response,  is  reported  to  the  outer 
loop  and  is  used  to  improve  the  estimate  of  the  impact  location. 


6.1  The  Inner  Loop  -  The  Smoot  her /Filter 

At  an  estimated  impact  location,  the  smoother/filter  reconstructs  a  force  history 
based  on  the  experimental  response.  This  technique  of  using  the  system  response  to 
find  the  system  properties  and  the  unknown  input  was  presented  by  Bryson  [11].  The 
technique  successfully  found  system  properties  of  a  helicopter  using  flight  test  data. 
The  following  development  closely  follows  Bryson’s  method. 
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The  reconstruction  of  the  force  history  begins  with  the  definition  of  the  figure-of- 
merit. 

J  =  ^  [zo  -  z(0)]^  So  [zo  -  z(0)]  +  7;Y^  /(^)Q/(”) 

^  (6.1) 

+  i  E  l/’■(n)R./(n)  +  -  [z„  -  z{NT)f  S,  -  z(NT)] 

^  .n=0  ^ 

where  Zo  and  z^t  are  user-supplied  guesses  of  initial  and  final  states.  /  is  the  re¬ 
constructed  force.  1/  is  the  error  between  the  experimental  data,  Sexp,  and  the  model 
response,  s.  The  error  contains  both  the  effects  of  measurement  noise  and  the  error 
due  to  incorrect  estimation  of  the  impact.  Q  is  the  diagonal  weighting  matrix  for  the 
input  force  and  R  for  the  states.  So  and  S /  are  weighting  matrices  for  the  initial  and 
final  conditions. 

The  values  Q  and  R  were  set  to  .05  and  5.0  x  10^®.  These  values  were  based  on 
Choi’s  [12]  work  on  the  beam  and  were  adjusted  empirically  for  the  composite  plate. 
Additional  information  on  these  values  can  be  found  in  Sage  and  Melsa  [13].  The 
system  is  at  rest  prior  to  the  impact  so  Zo  is  set  to  zero  and  So  is  set  to  a  very  high 
value.  Since  the  final  condition  is  unpredictable,  S/  is  set  to  zero. 

The  objective  is  to  minimize  J.  This  means  that  the  error,  u  ,  is  minimized,  the 
force  reconstruction,  /,  is  constrained  and  the  initial  conditions  are  constrained.  The 
force  is  unknown  but  its  reconstruction  is  constrained  by  penalizing  J  for  large  force 
values.  Similarly,  J  is  penalized  when  the  initial  condition  is  not  met. 

The  minimization  of  the  figure-of-merit  is  subject  to  the  system  equations  of 
Chapter  5,  Eq.  5.26  and  Eq.  5.27.  For  n  =  0,  . . . ,  {NT  —  1) 


z(n  +  1)  =  ^z(n)  -1-  r/(n) 

(6.2) 

s(n)  =  Cz(n) 

(6.3) 

where  NT  is  the  number  of  time  steps. 

To  satisfy  these  equations,  the  figure-of-merit  is  modified  to  accommodate  the 
constraints  using  Lagrange  multipliers  A 
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NT-1 

J  =  J  +  ^  A^(n  +  1)  [$z(n)  +  r/(n)  -  z(n  +  1)]  (6.4) 

n=0 

where  the  Lagrange  multiplier  was  introduced  at  an  extra  time  step,  n  =  NT.  For 
simplicity,  X{NT)  is  set  equal  to  zero. 

Minimizing  J  is  equivalent  to  minimizing  J  subject  to  the  constraint  equations. 
Taking  a  variation  on  the  modified  figure  of  merit  gives 

NT  NT-l 

w = E  +  E  («-5) 

Si  "/("J  ^ 


where  x  is  the  spatial  variable 


(6.6) 


At  a  stationary  point  of  J,  the  first  variation  vanishes  for  arbitrary  variation  of 
6z(n),  Sf(n)  and  6x,  and  so  the  coefficients  of  each  term  must  vanish.  The  first  two 
coefficients  are  solved  with  the  smoother/filter  while  the  third  coefficient  is  left  to  the 
outer  loop.  The  first  coefficient  is 


[zo-z(0)f So+A’’(l)$  forn  =  0 

dz(n)  I  — i/^(n)RC+A^(n  +  l)^-A'^(n)  forn  =  1,  . . . ,  NT 


The  second  coefficient  is 


= /^(n)Q+A^(n  +  l)r  for  n=  0,  ...,  (ATT  -  1) 

9f(n) 


(6.8) 


Setting  these  two  coefficients  equal  to  zero  and  incorporating  the  system  equations 
leads  to  the  smoothing  problem  for  a  given  x. 


z(n  +  1)  =  <&z(n)  +  r/(n) 

A(n)  =  $^A(n  +  l)-C^Ri/(n) 
f(n)  =  -Q-'r^A(n  +  l) 


(6.9) 
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The  boundary  conditions  are 


z(0)  =  Zo-So'^^A(l) 
A(iVr)  =  0 


(6.10) 


This  is  a  two-point  boundary  value  problem.  Half  of  the  boundary  conditions  are 
given  at  the  initial  point,  the  other  half  at  the  final  point.  The  boundary  value 
problem  can  be  solved  by  the  sweep  method.  The  final  boundary  condition  may  be 
swept  backward,  and  the  initial  condition  swept  forward. 

The  boundary  condition  at  n  =  NT  suggests  a  backward  sweep  solution  of  the 
form 


A(n)  =  S®(n)z(n)-A®(n)  (6-11) 

The  proper  form  of  S®  and  A^  need  to  be  found  to  solve  this  problem.  Using  the 
following  equation,  z(n  -f  1)  can  be  eliminated  from  Eq.  6.9. 


X(n) 

C^RC  -h  W„ 

w./ 

z(n) 

0 

1 

- 1 

_f(n) 

where 


rr 


A(71  -|-  1) 


C^R 

0 


(6.12) 


W„  W,f 
Wlf  Wff_ 


0  0 
0  Q 


rT 


S^(n-Hl)  [  ^  r  ] 


(6.13) 


Solving  the  second  equation  of  Eq.  6.12  for  /(n) 


/(n)  =  -WJ}  [W^/z(n)  -  TX^in  +  1)]  (6.14) 

and  substituting  the  result  into  the  first  equation,  results  in 


A(k)  =  [C’'RC  +  S'(n)]  z{n)  -  A'’(n)  -  C^Rv(n) 


(6.15) 
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where  the  time  downdate  relations  are  defined 


S5(n)  =  W,,-W,fWj}W^y 

A®(n)  =  -  W,/W7/r^]  A®(n  +  1) 


(6.16) 


This  time  downdate  process  involves  only  the  constraint  equations. 

Comparing  Eq.  6.11  and  Eq.  6.15,  the  following  measurement  downdate  relation 
can  be  defined 


S^(n)  =  S®(n)  +  C^RC 

A®(n)  =  X®(n)  +  C^'Rj^Cn) 


(6.17) 


This  step  refines  the  results  of  the  time  downdate  procedure  with  measurement  data. 

A  new  solution  method  is  presented  in  this  research  to  save  computation  time 
and  memory  requirements  associated  with  the  large  rank  one  matrix,  Bryson’s 
method  is  included  in  Appendix  B.  The  time  downdates  and  measurement  downdates 

for  and  are  combined  to  form  the  equations 

/ 

C®(n)  =  C®(n)$ 

Sf/n)  =  +  Sfy(n  +  1)  (6-18) 

Sff(n)  =  r^C^^RC®r  +  Q  +  Sf^(n+l) 


The  backward  sweep  starts  with  the  final  condition  and  the  last  recorded  data  point, 

Sexp(NT). 

C^(NT)  =  C 

AbCATT)  =  C®^(n)Rse.p(Arr) 

W^/(l)  =  $^C®^(n)RC^(n)r 
W//(l)  =  r^C®^(n)RC^(n)r 

The  sweep  proceeds  as  follows,  for  n  =  (iVT  —  1),  . .  • ,  1 
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KB(n)  =  Wj.}(NT-n)Wj^(NT-n) 
fain)  =  Wj}(NT-n)r^A^(n  +  l) 

\^{n)  =  ^^\^{n+.l)-W,f{NT-n)fB{n) 

\B(n)  =  A^(n)  +  C^Rs„p(n)  (6-20) 

C^(n)  =  C^(n)$ 

Sf^  =  +  Sf/n  +  1) 

sfj  =  r^c^^Rc^r  +  sf/n  + 1) 

Where  KB(n)  and  /b(w)  are  stored  for  use  in  the  forward  sweep. 

The  matrices  W^f  and  W//  require  a  new  updating  procedure  in  the  backward 
sweep.  Wii(n)  is  not  used.  W2/(n)  and  W//(n)  are  a  function  of  all  the  previous 
vectors.  For  i  =  1,  . . . ,  {NT  —  n) 

'  W2/(iVT  -  n  +  1)  =  W,f{NT-n  +  l)-^'^W,f{i)Wj}{i)Wjf{i)T 
Wff{NT-n  +  l)  =  WffiNT-n  +  l)-T^W,f{i)Wj}{i)Wjf{i)T 

The  S®  matrices  are  added. 

W,f{NT-n  +  l)  =  -W,f{NT-n+l)  +  Sff{n) 

WffiNT  -n  +  1)  =  Wff{NT  -n  +  l)  +  Sff{n) 

All  of  the  W  must  be  updated  for  the  next  step  in  the  backward  sweep. 

0,  ...,  (ATT-n-l) 

W2/(i)  =  ^W,f{l) 


(6.21) 


(6.22) 

For  i  = 

(6.23) 
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The  forward  sequence  is  unchanged  from  Bryson’s  method. 

For  n  =  1,  . . . ,  {NT  —  1) 

/(n)  =  /^(n)  -  K^(n)z(n) 

s(n)  =  Cz(n)  (6.24) 

z(n  +  1)  =  $z(n)  +  r/(n) 

/  is  the  reconstructed  force  history. 

Once  the  smoothing  problem  is  solved,  the  first  two  coefficients  of  Eq.  6.5  vanish, 
but  the  third  remains 

6J  =  =  jx<5x  (6.25) 

ox. 

The  third  coefficient  is 

Q  f  N-l  ■  Qp  N-l 

This  term  also  must  vanish  at  a  stationary  point  of  J.  This  is  done  in  the  outer  loop 
wh0re  the  gradient  information,  Ji,  is  used  to  improve  the  estimate  of  the  impact 
location  until  a  stationary  point  is  reached.  The  details  of  the  outer  loop  that  solves 
this  non-linear  location  identification  problem  is  discussed  in  Section  6.3. 

Since  A,  i/,  /  and  z  are  all  calculated  for  the  other  two  coefficients,  this  third  coeffi¬ 
cient  can  be  readily  obtained.  This  is  one  of  the  great  benefits  of  this  smoother/filter; 
it  provides  the  information  necessary  to  update  the  estimate  of  the  location. 


6.2  Advantages  of  New  Filter 

The  new  backward  sweep  method  offers  significant  computational  savings  for  sys¬ 
tems  with  a  large  number  of  degrees  of  freedom.  The  savings  in  required  memory  is 
described  below.  Similar  savings  are  seen  with  the  multiplication  requirements. 

The  elenicnt  of  Bryson’s  backward  filter  that  consumes  most  of  the  memory  stor¬ 
age  is  the  matrix. 
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In  the  present  application  there  is  a  single  force  and  N^of  degrees  of  freedom. 
Where  N^of  =  M  x  N  (the  number  of  x  and  y  frequency  components).  The 
matrix  then  has  [Ndof  x  I^dof]  elements. 

In  the  new  formulation  the  large  rank-one  matrix  is  replaced  with  a  series  of 
vectors.  The  element  of  the  filter  that  uses  the  most  memory  storage  is  the  matrix 
multiplication  in  Sf^ 

$^w,/(0W7/(i)wf/i)r  (6.27) 

The  size  of  the  block  diagonal  matrix  $  is  [Ndof  x  Ndof]  and  F  is  [Ndof  x  1].  The 
multiplication  results  in  a  Sff  matrix  of  size  [Ndof  x  1].  One  of  these  vectors  is  needed 
for  each  time  step.  The  storage  requirements  are  for  [Ndof  x  NT]  matrix. 

In  summary,  the  largest  contributer  to  memory  requirements  for  the  new  method 
is  proportional  to  Ndof  x  NT,  as  compared  to  N^of  for  Bach’s  method. 

In  a  problem  with  many  degrees  of  freedom,  as  in  the  impact  on  a  composite 
plate,  the  new  method  offers  significant  computational  savings.  The  exact  memory 
and  computational  time  requirements  of  the  new  method  and  Bryson’s  method  are 
compared  in  Table  1.  With  the  computational  savings,  it  is  possible  to  attain  near 

real-time  application  of  the  identification  system  to  the  composite  plate. 

/ 

/ 

6.3  The  Outer  Loop  -  The  Estimate  of  the  Impact 
Location 

The  outer  loop  of  the  inverse  problem  solver  is  the  search  for  the  impact  loca¬ 
tion.  The  search  finds  the  minimum  of  the  figure-of-merit  with  as  few  calculations  as 
possible  of  the  filter/smoother. 

The  estimate  of  the  impact  location  is  updated  using  the  gradient  of  the  figure 
of  merit  J  calculated  in  the  smoother/filter.  To  differentiate  between  the  local  and 
global  minima  of  this  nonlinear  problem,  the  search  starts  from  multiple  initial  guesses 
and  continues  with  a  series  of  line  searches.  Among  the  minima  found,  the  one  with 
the  lowest  value  of  the  figure-of-merit  is  taken  as  the  global  minimum. 

The  initial  guesses  are  distributed  throughout  the  area  that  is  being  monitored  for 
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an  impact.  From  each  of  these  guesses  a  line  search  will  start  using  the  quasi-Newton 
procedure 


Xe(i  +  1)  =  Xe(i)  -  ji(i)  (6.28) 

where  Xg  is  the  estimate  of  the  impact  location,  jii(i)  is  the  Hessian  matrix,  and  a 
is  the  step  size  scaling  variable,  i  is  the  counter  for  each  line  search. 

The  first  line  search  starts  in  the  direction  of  steepest  descent,  =  1-  From 

the  minimum  of  this  line  search  another  line  search  starts.  The  Hessian  matrix  is 
estimated  numerically  from  successive  values  of  the  gradient  vector  Jj,  using  a  rank- 
two  update  procedure.  See  Appendix  C.  The  step  size  is  determined  by  a.  The  series 
of  searches  continues  until  a  minimum  of  the  figure  of  merit  is  reached,  when  Ji  =  0. 

Within  the  line  search,  multiple  evaluations  of  the  smoother/filter  are  performed 
to  find  the  minimum  of  the  figure-of-merit  along  the  line.  The  line  search  starts  by 
taking  the  maximum  step  size  in  the  direction  determined  from  Eq.  6.28.  The  next 
step  will  be  along  the  line  but  half  the  size  of  the  first,  either  forward  or  backward. 
The  direction  is  determined  by 


jx(0)  •  jx(j) 

Jx 


(6.29) 


where 


ji(i)  =  gradient  at  current  point  of  line  search 
Jj.(0)  =  gradient  at  initial  point  of  line  search  (6.30) 

I  ji|  =  the  magnitude  of  the  vector 


and  j  is  the  counter  within  the  line  search. 

If  K  >  0,  then  the  search  has  not  reached  the  minimum,  and  the  search  continues 
forward.  If  «  <  0,  then  the  minimum  has  been  passed,  the  gradient  has  changed 
sign,  and  the  search  takes  a  step  back.  The  minimum  of  the  line  search  will  then  be 
reached  after  several  evaluations  of  the  smoother /filter,  see  Figure  13. 
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The  maximum  error,  eu,  and  range,  L,  of  the  line  search  can  be  calculated  from 


e/s 


I 

2nu 


(6.31) 


L 


/  {  1  + 


2n(.-l  _  1\ 

j 


where 


I  =  maximum  step  size  of  line  search 
nu  —  number  of  smoother/filter  evaluations  per  line  search 


(6.32) 


The  objective  of  the  search  technique  is  to  minimize  the  number  of  times  the 
smoother/filter  is  used  to  update  the  estimated  location  while  searching  for  the  global 
minimum  in  the  figure  of  merit.  Several  tests  were  used  to  control  the  start  and 
continuation  of  each  line  search.  The  following  tests  were  designed  to  avoid  searching 
in  an  area  far  from  the  impact  and  to  finish  the  search  close  to  a  minimum.  Because 
experimental  data  is  used,  it  is  not  practical  to  look  for  an  exact  zero  in  the  gradients 
of  the  figure-of-merit  contour.  If  the  gradients  are  less  than  a  preset  tolerance,  they 
are  considered  zero. 

•  If  the  gradient  of  the  figure-of-merit,  Jj,  is  near  zero  at  the  initial  point  of  a  line 
search,  the  search  is  cancelled.  Either  the  estimated  location  is  very  far  from 
the  impact  or  is  already  at  a  minimum.  The  value  of  J  is  stored  as  a  possible 
global  minimum. 

•  If  the  value  of  J  at  the  end  of  the  line  search  has  not  decreased  by  10%  or 
more  from  the  starting  point  of  the  search,  then  the  search  is  cancelled.  This  is 
an  additional  safeguard  against  searching  in  flat  regions  of  the  figure  of  merit 
contour. 

•  If  the  line  search  is  perpendicular  to  the  direction  of  steepest  descent,  the  search 
is  cancelled.  This  allows  for  an  early  exit  from  the  line  search  if  the  minimum 
is  reached  before  the  maximum  number  of  evaluations. 
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Figure  13:  The  line  search.  The  line  search  continues  through  a  series  of 
branch  points.  The  gradient  of  the  figure-of-merit  is  used  to  guide 
the  search.  The  search  is  ended  when  the  gradiant  of  J  is  less  then 
the  preset  tolerance. 
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Table  1:  The  computational  advantages  of  the  new  backward  filter  for  sys¬ 
tems  with  a  large  number  of  degrees  of  freedom.  Ndof  is  the  number 
of  degrees  of  freedom  in  the  system  and  NT  is  the  number  of  time 
steps. 


Number  of  multiplies 
Elements  to  store 


Present  Method 
8  X  Ndof  X  Np 
10  X  Ndof  X  NT 


Previous  Method 


3  X  N^^f  X  NT 


3xAr,V 


Chapter  7 

Computer  Code  -  IDIMPACT 


With  the  theory  of  the  identification  system  developed  in  the  previous  chapters, 
the  system  was  implemented  with  the  computer  code  IDIMPACT.  As  indicated  in 
Figure  14,  IDIMPACT  uses  the  digitized  plate  reponse  and  application  specific  pa¬ 
rameters  to  identify  the  impact  and  report  the  results  in  a  graphic  interface.  The 
code  provides  near  real-time  impact  identification. 

The  next  chapter  presents  the  application  of  IDIMPACT  to  a  composite  plate. 
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Figure  14:  The  input  and  output  of  IDIMPACT. 


Chapter  8 

Application  of  IDIMPACT 


With  the  theoretical  development  and  computer  implementation  of  the  identifica¬ 
tion  system,  the  next  step  is  designing  the  system  for  a  specific  application.  Appli¬ 
cation  parameters  are  specified  for.each  component  of  the  identification  system.  The 
system  was  applied  to  the  TSOO  composite  plate  described  in  the  Chapter  5. 

The  objective  in  designing  the  system  is  to  keep  the  total  solution  time  of  the 
force  identification  to  a  minimum.  The  solution  time  is  consumed  mostly  by  the 
smbother/filter  in  the  inverse  problem  solver.  To  keep  the  solution  time  down,  the 
solution  time  of  the  smoother/filter  and  the  number  of  times  the  filter  is  called  must 
be  kept  to  a  minimum.  The  speed  of  the  smoother  depends  on  the  forward  model. 
The  number  of  times  the  smoother  is  called  depends  on  the  outer  loop  of  the  inverse 
problem  solver  -  the  search  for  the  impact  location. 

An  example  impact  was  used  to  demonstrate  all  of  the  components  of  the  identi¬ 
fication  system. 


8.1  Experimental  Setup 

The  experimental  setup  includes  the  composite  plate,  sensor  array,  instrumented 
hammer,  data  acquisition  system  and  computer  controller  and  processor.  The  system 
was  designed  so  that  the  data  could  be  acquired  and  processed  using  one  computer. 
The  following  hardware  and  software  items  were  used: 
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•  r800[45/90/  -  45/0/0/45/0/0/  -  45/0]s  carbon  composite  plate  described  in 
Chapter  5.  In  this  application  the  plate  had  free  boundaries  with  support  only 
in  the  corners.  The  example  impact  was  at  {x,y)  =  (12.5,5.0). 

•  Piezo  Kinetics  Incorporated,  PKI-400  .25  inch  diameter  disk,  piezoelectric  sen¬ 
sors  discussed  in  Chapter  4.  The  response  of  the  plate  to  the  impact  was  mea¬ 
sured  using  a  thirteen  sensor  array  of  piezoelectrics  bonded  to  the  composite 
plate.  An  area  of  20  inches  by  20  inches  was  covered.  The  geometry  is  shown 
in  Figure  15.  Each  sensor  was  grounded  to  the  composite  plate. 

•  PCB  Piezotronics  instrumented  hammer,  208  A04  force  transducer  and  480A 
power  unit.  The  output  of  the  load  cell  was  recorded  for  comparison  to  the 
reconstructed  force. 

•  National  Semiconductor  LF  412  operational  amplifier.  The  amplifiers  were  used 
for  the  sensors  and  the  instrumented  hammer. 

•  Kiethley  Metrabyte  DAS  1800  analog  to  digital  signal  converters.  Each  board 
has  eight  channels  with  a  maximum  conversion  rate  of  333Msamples/sec. 

/ 

•  NCA  133MHz  pentium  computer. 

•  Keithley  Metrabyte  VTX  controller  for  the  DAS1800. 

•  Microsoft  Windows  95. 

•  Microsoft  C-l— I-  and  Visual  Basic. 

The  data  acquisition  system  was  controlled  by  the  VTX  controller  operating  in 
Visual  Basic.  It  was  triggered  by  the  voltage  output  of  the  instrumented  hammer. 
When  the  signal  from  the  hammer  crosses  a  predetermined  threshold,  a  TTL  trigger 
signal  is  sent  to  the  data  acquisitions  boards.  The  voltage  was  amplified  to  reduce 
noise  and  to  provide  a  zero  offset.  The  hammer  amplifier  is  described  in  Figure  17. 

Once  the  data  acquisition  is  triggered  the  voltage  generated  by  the  sensors  is 
digitized  at  22.3  KSamples/sec.  150  data  points  from  each  channel  are  recorded  -  10 
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pre-trigger  and  140  post  trigger.  The  piezoelectric  sensors  are  very  sensitive  to  strain 
and  unamplified  easily  produce  output  over  ten  volts.  Amplifiers  provided  a  zero 
offset  for  the  sensor  voltage  output.  The  sensor  amplifier  is  described  in  Figure  16. 

Once  the  data  is  digitized,  processing  can  begin.  IDIMPACT,  described  in  Chap¬ 
ter  7,  operates  in  the  32  bit  environment  of  Windows  95. 

8.2  Sensor  and  Preprocessor 

The  gain  of  the  sensors  was  determined  by  comparing  the  response  of  the  forward 
niodel  with  the  measured  response,  ge  =  4.0  x  10“®(l/voZfs). 

The  preprocessor  divided  the  full  20  inch  by  20  inch  sensor  array  into  five  overlap¬ 
ping  regions  each  with  five  sensors.  (Five  sensors  were  used  NS  =  5.)  After  scanning 
the  sensors  to  determine  which  had  the  highest  initial  energy,  one  of  the  regions  is  cho¬ 
sen  to  be  used  by  the  inverse  problem  solver.  The  geometry  of  the  isolation  regions, 
the  definition  of  the  global  and  local  coordinate  systems,  and  the  sensor  numbering 
system  are  shown  in  Figure  18. 

With  the  example  impact  at  {x,  y)  =  (12.5, 5.0),  the  region  of  five  sensors  with  the 
highest  initial  energy  was  the  the  lower-right  region,  region  four.  The  sensor  readings 
at  these  five  sensors  is  shown  in  Figure  19.  In  the  local  coordinates  of  isolation  region 
four,  the  impact  is  at  {x',y')  =  (2.5, 5.0) 

The  different  energy  content  in  the  initial  time  window  of  the  sensors  at  different 
distances  from  the  example  impact  is  shown  in  Figure  20.  With  the  example  impact, 
sensor  4  located  at  {x,y)  =  (10,0)  is  5.59  inches  from  the  impact  while  sensor  5 
located  at  (x,y)  =  (20,0)  is  13.5  inches  from  the  impact.  The  time  window,  selected 
so  the  sensor  farther  away  would  have  a  very  small  energy  content  compared  to  the 
close  sensor,  was  26.9  x  10~®  seconds.  This  corresponds  to  the  first  seven  data  points 
in  the  sensor  data.  Ne  =  7. 
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Figure  15:  The  composite  plate  with  the  dimensions  and  geometry  of  the 
thirteen  piezoelectric  sensor  array.  The  global  coordinate  system 
is  shown.  The  sensors  are  on  the  back  of  the  plate.  The  plate 
with  free  boundary  conditions  is  shown. 
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Force  Transducer 
PCB  208  A04 


Figure  16:  Wire  diagram  of  hammer  amplifier. 
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Amplifier 


Figure  17:  Wire  diagram  of  sensor  amplifier.  No  amplification  is  applied  but 
the  operational  amplifier  is  used  to  provide  a  zero  offset  to  the 
signal. 
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Figure  18:  The  sensor  array  and  the  five  overlapping  isolation  regions,  Rl- 
R5.  The  global  coordinates  (a;,  y)  and  the  local  isolation  region 
coordinates  {x',  y')  are  shown. 
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Figure  19:  The  response  of  the  five  sensors  in  isolation  region  four.  The 
impact  was  reported  at  {x',  y')  =  (2.7, 5.2).  The  actual  impact 
was  at  (x',  y')  =  (2.5, 5.0). 
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Figure  20:  Response  of  the  sensor  at  (x',  y')  =  (0.0, 0.0)  (solid  line)  and  at 
{x',y')  =  (10.0,0.0)  (dashed  line)  to  the  example  impact.  The 
energy  window  is  used  to  determine  the  sensors  closest  to  the 
impact. 
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8.3  Forward  Model 

The  forward  model  must  accurately  estimate  the  plate  response  while  using  as 
few  time  steps  and  frequency  components  as  possible.  The  time  steps  and  frequency 
components  determine  the  solution  time  of  the  inverse  problem  solver. 

The  periodic  solution  of  the  Fourier  transform  can  lead  to  an  error,  aliasing, 
when  the  solution  overlaps.  It  is  necessary  to  have  a  large  spatial  domain  to  prevent 
overlapping  as  the  wave  propagates.  The  spatial  domain  was  set  to  X  =  F  =  90 
inches.  To  maintain  the  spatial  resolution,  many  frequency  components  were  needed, 
M  =  N  =  180. 

The  time  duration  of  the  measured  data  has  to  be  long  enough  to  include  the  full 
impact  force  history.  The  duration  of  the  impact  was  approximately  one  millisecond. 
Two  milliseconds  of  data  was  recorded  to  insure  all  of  the  force  was  included.  The 
number  of  time  steps,  NT,  was  45,  the  time  step,  Ts,  was  4.48  x  10"®  seconds.  This 
results  in  a  total  time  duration  of  Ts  x  {NT  —  1)  =  1-97  milliseconds. 

For  this  application,  it  was  discovered  that  the  frequency  components  beyond  the 
first  36  were  essentially  zero.  Instead  of  using  all  the  components  in  the  solution,  the 

144  i^ear  zero  components  were  treated  as  zero  padding.  A  great  deal  of  computation 

! 

tiifie  is  saved  using  only  the  non-zero  components. 


8.4  Inverse  Problem  Solver 

The  outer  loop  of  the  inverse  problem  solver,  the  search  for  the  impact  location, 
must  be  designed  so  that  the  impact  is  located  with  as  few  evaluations  of  the  inner 
loop,  the  smoother/filter,  as  possible  while  still  maintaining  the  robustness  of  the 
system. 

•  The  weighting  values  for  the  smoother/filter  were  determined  in  Chapter  6.  The 
values  Q  and  R  were  set  to  .05  and  5.0  x  10^®. 

To  help  design  the  search  technique,  an  exhaustive  search  was  performed  for  the 
example  impact.  The  smoother/filter  calculated  the  figure  of  merit,  J,  on  a  grid  with 
1/2  inch  spScing  in  the  isolation  region.  This  creates  a  figure-of-merit  contour.  The 
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minimum  of  this  contour  is  the  identification  system’s  estimate  of  the  actual  impact 
location.  The  typical  contour  has  one  global  minimum  and  several  local  minima,  see 
Figure  21. 

There  are  three  components  of  the  search  technique  to  be  designed  for  this  contour, 
the  number  and  spacing  of  the  initial  guesses,  the  accuracy  and  range  of  the  line 
search,  and  the  tests  to  control  the  start  and  continuation  of  each  line  search. 

There  has  to  be  enough  initial  guesses  distributed  over  the  search  region  to  assure 
that  the  global  minimum  can  be  found.  It  was  observed  from  experiment  that  with  an 
impact  near  the  edge  of  the  sensor  array,  it  was  difficult  to  find  the  contour  leading 
to  the  mininnim  of  the  figure-of-merit.  With  an  impact  near  the  middle,  it  was 
relatively  easy  to  find  the  correct  contour.  To  assure  that  the  global  minimum  could 
be  found,  the  initial  guesses  were  distributed  near  the  search  region  boundary.  The 
initial  guesses  are  shown  in  Figure  18. 

To  maintain  a  robust  search,  there  should  be  an  overlap  in  the  range  of  each 
search.  The  maximum  step  size,  I,  of  the  first  line  search  was  chosen  as  1.6  inches 
with  3  evaluations  of  the  filter/smoother,  riis  =  3.  The  successive  line  searches  were 
limited  to  a  step  size  of  1.0  inches  and  3  function  evaluations.  Using  Eq.  6.31,  the 
range  for  the  first  line  search  is  2.80  inches  and  for  the  successive  searches  is  1.5 
inches.  Since  three  successive  line  searches  are  allowed,  the  maximum  total  range 
of  each  search  is  5.8  inches.  The  maximum  error  for  the  minimum  of  the  first  and 
successive  line  searches  is  .27  and  .25  inches.  The  accuracy  and  overlapping  searches 
provided  good  results  for  the  application  to  the  T800  composite  plate. 

The  preset  tolerance  for  the  gradiant  needed  to  start  a  search  was  Je  =  40.  The 
tolerance  needed  to  approximate  a  minimum  of  the  figure-of-merit  was  Jg  =  15.  If 
Je  <  40,  the  search  was  not  started  firom  the  initial  guess.  If  Jg  <  15,  the  search  was 
stopped  because  a  minimum  was  found. 

The  new  formulation  of  the  smoother /filter  offers  a  great  computational  advantage 
for  this  application.  With  the  number  of  non-zero  frequency  components  and  time 
steps  in  the  forward  model,  there  is  a  91%  saving  in  the  computational  time  and 
an  88%  savings  in  the  memory  requirements  for  the  new  filter,  see  Table  2.  The 
computational  savings  were  calculated  using  the  relations  in  Table  1.  (Eight  bytes 
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Figure  21:  The  figure  of  merit  contour.  The  smoother /filter  evaluated  the 
figure-of-merit  on  a  1/2  inch  grid  in  the  isolation  region.  The 
actual  impact  was  at  =  (2.5, 5.0) 
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are  stored  per  double  on  the  personal  computer.) 

8.5  Identification  of  the  Example  Impact 

With  all  of  the  application  parameters  determined,  IDIMPACT  was  run  to  identify 
the  example  impact.  Prom  two  starting  points,  the  impact  was  found,  see  Figure  22. 
The  two  line  search  stopped  at  slightly  different  points  because  the  searches  are 
limited  to  three  consecutive  line  searches.  The  searches  starting  from  the  upper  right 
and  lower  right  corner  were  cancelled  because  the  value  of  the  figure  of  merit  at  the 
end  of  the  line  search  had  decreased  by  less  than  10%  from  the  starting  point.  The 
searches  were  not  started  from  the  other  initial  guesses  because  the  gradiant  was 
below  the  preset  tolerance,  Jg  <  40.  The  impact  was  reported  .28  inches  from  the 
actual  location  and  the  total  impact  energy  reported  was  within  12%  of  the  actual. 

A  total  of  32  evaluations  of  the  smoother /filter  was  needed  to  complete  the  search. 
On  the  PC,  each  function  evaluation  took  8  seconds  for  a  total  solution  time  for  the 
identification  system  of  4  minutes  and  16  seconds.  The  performance  of  the  identifi¬ 
cation  system  allows  near-real  time  application  of  IDIMPACT . 

''The  reconstructed  force  is  compared  to  the  experimentally  recorded  force  in  Fig¬ 


ure  23. 
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Figure  22:  Search  using  multiple  starting  points.  The  search  starting  from 
(x',  y')  =  (1.6, 5.0)  successfully  found  the  impact  location.  The  re¬ 
ported  impact  location  was  at  (x',  y')  =  (2.7, 5.2)  while  the  actual 
location  was  {x',  y')  =  (2.5, 5.0). 
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Figure  23:  The  reconstructed  force  compared  to  the  example  impact. 
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Table  2:  The  new  filter  applied  to  the  plate  impact  identification  problem. 

To  model  the  plate  response,  36  frequency  components  in  the  x  and 
y  direction  were  used  with  45  time  steps. 


Present  Method 

Previous  Method 

Percent  Savings 

Number  of  multiplies 

21  Million 

227  Million 

91 

Memory  Required 

4.7  MBytes 

40.3  MBytes 

88 

Chapter  9 

Identification  System  Accuracy 
and  Robustness 


In  the  previous  chapter  the  identification  system  was  designed  and  applied  to  the 
T800  composite  plate.  The  system  successfully  identified  the  example  impact.  In  this 
chapter  the  identification  system  will  be  tested  with  numerous  impacts  to  quantify 
the  error  in  its  estimate  of  the  impact  location  and  the  force  reconstruction.  The 
system  will  be  tested  on  the  plate  with  clamped  and  free  boundary  conditions.  Also, 
to  further  test  the  robustness  of  the  system,  a  noisy  environment  will  be  simulated. 


9.1  The  Test  Impact  Array 

An  array  of  impacts  was  applied  to  the  composite  plate.  The  distribution  of  40 
impacts  is  shown  in  Figure  24.  The  impacts  were  distributed  over  a  quarter  of  the 
plate,  isolation  region  4.  Because  the  plate  is  symmetric,  the  results  from  the  region 
represent  whole  sensor  array. 

A  direct  hit  over  the  sensor  was  not  included  in  the  impact  distribution  due  to 
a  limitation  of  the  data  acquisition  system.  With  an  impact  very  near  a  sensor,  the 
sensor  output  will  exceed  10  volts.  This  saturates  the  sensor  channel  on  the  analog  to 
digital  converter.  The  other  sensor  channels  are  also  affected  and  the  digitized  data 
is  not  accurate.  With  the  saturated  data  and  the  inaccurate  data,  it  is  not  possible 
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for  the  identification  system  to  identify  the  impact.  This  limitation  can  be  overcome 
with  improvement  in  the  acquisition  system;  it  is  not  a  theoretical  limitation  of  the 
identification  system. 


9.2  Free  Plate  Boundaries 

The  impact  test  array  was  applied  to  the  composite  plate  with  the  free  boundary 
conditions.  The  reported  location  results  are  shown  in  Figure  25,  where  the  actual 
impact  is  at  the  center  of  the  circle  and  the  reported  impact  is  in  error  in  the  x  or 
p  direction.  The  reported  impact  location  was,  on  the  average,  0.39  inches  from  the 
actual  impact  location.  71.1%  of  the  impacts  were  within  0.5  inches  of  the  actual. 

The  error  in  the  energy  of  the  reconstructed  force  is  defined  as 


I]  fhpW 

n=0 


where  Brecon  and  Egxp  are  the  energy  of  the  reconstructed  force  and  the  experimentally 
recorded  force,  fexp  is  the  recorded  force  from  the  instrumented  hammer.  The  energy 
of  the  reconstructed  force  differed  by  an  average  of  14.2%  from  the  recorded  impact. 
The  average  error  in  the  peak  magnitude  was  3.1%.  The  reconstruction  of  the  example 
impact  was  shown  in  the  previous  chapter.  Figure  23. 

Noise  was  added  to  the  digitized  data  obtained  from  the  impact  test  array.  This 
simulated  noise  was  calculated  using 


y noise  =  Veip  +  0.5  X  Raudu 

^noise  ~  9e  ^  ^exp 


(9.2) 
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where  the  Randn  function  randomly  generates  a  normal  distribution  of  numbers  with 
a  mean  of  0.0  and  a  variance  of  1.0.  The  noise  generated  is  approximately  10%  of 
the  magnitude  of  the  original  experimental  data.  The  plate  response  to  the  example 
impact  with  added  noise  is  shown  in  Figure  26. 

The  reported  location  results  are  shown  in  Figure  27.  The  average  impact  reported 
for  the  noisy  case  was  .41  inches  from  the  actual.  73.7%  of  the  impact  were  located 
within  0.5  inches  of  the  actual.  The  example  impact  was  reported  at  {x',y')  = 
(3.0, 5.3). 

The  energy  of  the  reconstructed  force  history  was  on  the  average  within  11.8% 
of  the  recorded  impact  history.  The  average  error  in  the  peak  magnitude  was  2.0%. 
The  reconstruction  of  the  impact  at  the  example  location  is  shown  in  Figure  28. 

9.3  Clamped  Plate  Boundaries 

The  impact  array  was  also  used  to  test  the  identification  system  on  the  plate 
with  clamped  boundaries.  The  reported  location  results  are  shown  in  Figure  30.  The 
average  error  of  the  reported  impact  was  .47  inches.  66.7%  were  found  with  0.5  inches. 
The  example  impact  was  reported  at  {x',y')  =  (2.7,5.!). 

The  example  impact  is  reconstructed  in  Figure  31.  With  the  forty  impacts  tested, 
the  average  energy  of  the  reconstructed  force  was  within  13.4%  percent  of  the  ex¬ 
perimentally  recorded  force.  The  magnitude  of  the  force  was  reconstructed  within 
8.25%. 

Noise  was  also  added  to  the  data  recorded  with  the  clamped  boundary  conditions. 
The  results  are  in  Figure  33.  The  average  error  of  the  reported  impact  was  .49 
inches.  60.5%  were  found  with  0.5  inches.  The  example  impact  was  reported  at 
(x',y')  =  (2.7,5.!). 

The  example  impact  is  reconstructed  in  Figure  31.  The  average  energy  of  the 
reconstructed  force  was  within  16.3%  percent  of  the  experimentally  recorded  force. 
The  magnitude  of  the  force  was  reconstructed  within  9.5%. 
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maximum  error  is  1.28  inches 


Figure  25:  Results  from  ten  inch  sensor  spacing  with  free  plate  boundaries. 
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Milliseconds  Milliseconds 


Figure  26:  The  response  of  the  example  impact  with  noise  added  to  the  sen¬ 
sor  response.  The  plate  has  free  boundary  conditions. 
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maximum  error  is  1.12  inches 


Figure  27:  The  accuracy  of  the  impact  identification  system  with  the  noisy 
sensor  measurement  and  the  free  boundary  conditions. 
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Milliseconds 


Figure  28:  The  reconstructed  force  with  noise  in  the  recorded  sensor  data. 

The  example  impact  at  {x',y')  =  (2.5, 5.0)  was  reported  at 
(3.0, 5.3).  The  plate  has  free  boundary  conditions. 
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Figure  29:  The  response  of  the  five  sensors  to  an  impact  at  (x',  y')  = 
(2.5, 5.0).  The  boundaries  of  the  plate  are  clamped. 
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Milliseconds 


Figure  31:  The  typical  reconstructed  force  history  for  the  clamped  boundary. 

This  is  the  impact  reconstruction  for  the  impact  at  {x',y')  = 
(2.5, 5.0).  The  reported  impact  was  at  {x',y')  =  (2.7, 5.1) 
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Milliseconds  Milliseconds 


Figure  32:  The  response  of  the  five  sensors  to  an  impact  at  (x',  y')  = 
(2.5, 5.0).  The  boundaries  of  the  plate  are  clamped  and  noise 
is  added  to  the  sensor  data. 
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maximum  error  is  1.06  inches 


Figure  33:  Result  from  sensor  array  with  clamped  plate  boundaries  and  noise 
added  to  sensor  data. 
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Figure  34:  The  typical  reconstructed  force  history  for  the  clamped  boundary 
with  noise  added  to  the  sensor  data.  This  is  the  impact  recon¬ 
struction  for  the  impact  at  {x',y')  =  (2.5, 5.0).  The  reported 
impact  was  at  {x\y')  =  (2.7, 5.1) 
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9.4  Discussion  of  Results 

The  identification  proved  to  be  a  very  reliable  and  robust  system  for  locating  and 
reconstructing  the  force  time  history  of  an  unknown  impact.  The  system  had  very 
accurate  and  consistent  performance  in  the  four  cases  tested:  the  composite  plate 
with  free  or  clamped  boundary  conditions  and  with  clean  or  noisy  data.  The  average 
reported  location  differed  by  .39  inches  for  the  free  boundary  case  to  .49  inches  for 
the  clamped  plate  with  noise  added  to  the  sensor  data.  The  best  match  between 
the  modal  response  and  the  experimental  response  was  achieved  with  the  free  plate 
boundaries. 

The  most  significant  error  in  the  reported  location  occured  when  the  impact  was 
close  to  a  sensor  on  the  edge  of  the  sensor  array.  For  example,  the  largest  error  of 
the  clamped  boundary  case  was  1.06  inches.  This  occured  when  the  plate  was  hit  at 
{x',y')  =  (10,1.25).  The  nearby  sensor,  sensor  5  at  {x',y')  =  (10.0,0.0),  measures 
strain  of  significantly  greater  magnitude  then  the  other  sensors  in  the  isolation  region. 
In  the  search  for  the  impact  location,  this  sensor  has  a  big  effect  on  the  value  of  the 
figure-of-merit.  When  the  model  response  closely  matches  this  sensor  measurement, 
the:  figure-of-merit  has  a  value  near  its  minimum.  This  occurs  when  the  estimated 
location  is  the  correct  radius  from  the  close  sensor.  The  search  will  report  the  location 
at  this  radius  which  may  be  a  significant  distance  from  the  actual  impact. 

A  possible  way  to  correct  this  problem  is  to  scale  each  value  of  the  diagonal 
weighting  matrix  R  so  that  each  sensor  has  approximately  the  same  influence  on  the 
figure-of-merit.  The  values  of  R  could  be  scaled  by  comparing  the  magnitude  of  the 
data  of  the  five  sensor  in  the  isolation  region. 

The  identification  system  successfully  identified  the  impact  with  noise  added  to  the 
sensor  data.  The  noise  effected  the  accuracy  of  the  estimated  location  slightly.  The 
average  error  in  the  estimated  location  was  increased  by  5%  to  10%.  The  accuracy 
of  the  force  reconstruction  was  not  greatly  effected. 

The  reconstruction  provided  a  very  good  estimate  of  the  energy  of  the  impact  force. 
The  greatest  difference  between  the  recorded  and  the  reconstructed  force  occurred 
near  the  end  of  the  time  history.  This  difference  is  due  to  increasing  error  between 
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Chapter  10 
Conclusion 


The  impact  identification  system  proved  to  be  a  very  robust  and  accurate  means 
of  identifying  the  impact  on  a  composite  plate.  The  system  was  tested  with  many 
impacts  distributed  over  the  composite  plate  with  both  clean  and  noisy  sensor  data 
and  with  free  and  clamped  boundaries.  The  system  demonstrated  the  effectiveness 
of  the  identification  system  based  on  the  smoother/filter. 

There  were  key  developments  that  made  the  identification  system  applicable  to 
structures. 

•  The  derivation  of  the  governing  equation  for  a  composite  plate  in  state  space 
form.  This  two  dimensional  solution  was  extended  from  a  one  dimensional  beam 
solution. 

•  The  reformulation  of  the  filter /smoother  to  handle  systems  with  a  large  number 
of  degrees  of  freedom.  The  reformulation  saved  considerable  computation  time 
in  this  structures  application. 

•  The  use  of  small  circular  piezo-electric  sensors  to  measure  a  non-directional 
strain  property  on  the  surface  of  the  composite  plate.  These  sensors  provided 
a  simple  and  reliable  means  of  measuring  the  plate  response. 

Future  work  could  include  the  application  of  the  identification  system  to  more 
complicated  problems.  The  forward  model  can  be  extended  to  include  more  general 
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structures.  Finite  element  models  can  also  be  used.  The  inverse  problem  solver  based 
on  the  smoother/filter,  is  a  very  general  solution  technique  and  can  be  applied  to 
other  new  structures  problems. 
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Appendix  A 

The  Observation  Matrix 


The  calculation  of  the  plate  deflection  and  strain,  presented  in  Chapter  5,  uses 
a  simplified  submatrix  C^,,  of  the  observation  matrix,  C.  This  simplified  matrix  is 
attained  by  taking  advantage  of  the  symmetry  of  the  impact  problem. 

The  equation  used  to  calculate  the  plate  deflection,  w,  begins  with  the  two- 
dimensional  Fourier  Transform  between  the  space  and  frequency  domain. 

iV-lM-l 

n=0  m=0 

where  Wmn  are  the  frequency  components  of  the  plate  displacement.  Km  = 
and  Kn  =  M  and  N  are  the  number  of  frequency  components  in  the  x  and  y 
directions.  X  and  Y  are  the  extent  of  the  spatial  domain.  This  equation  can  be 
written  in  matrix  form 

•w{xi,yi,t)  =  C^u>  (A.2) 

where  a>  is  a  [A  •  M  x  1]  vector  with  elements 

u>(n  +  Nm)  =  Umn  (A.3) 

and  Cu,  is  a  [1  X  M  •  A]  vector  for  each  sensor  with  elements 

C„(n  -1-  Am)  =  (A.4) 

where  m  =  0,  . . . ,  (M  —  1)  and  ra  =  0,  . . . ,  (A  —  1)  . 
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To  better  visualize  the  symmetry  of  the  problem,  the  matrix  is  temporarily 
presented  as  a  [M  x  N]  matix. 


C^(nm)  =  (A.5) 

The  size  of  this  matrix  can  be  reduced  to  approximately  one  quarter  of  its  original 
size  with  the  symmetric  properties  of  the  exponents 


(A.6) 


The  reduction  in  the  size  of  the  can  be  visualized  by  dividing  the  matrix  into 
four  sections  about  the  m  =  y  and  n  =  y  element.  The  matrix  can  then  be  folded 
about  the  x  and  p  divisions.  The  resulting  elements,  with  the  exception  of  the  edge 
elements,  are  the  addition  of  the  four  sections. 

The  matrix  can  be  further  simplified.  The  displacement  of  the  plate,  w,  is  a  real 
and  symmetric  function  about  the  x  and  y  axes.  The  Fourier  Transform  of  a  real  and 
even  faction  in  the  space  domian  results  in  a  real  and  even  faction  in  the  frequency 
domain.  This  means  the  components  of  (Cu,)mn  are  real  and  symmetric  about  the 
n  A  A  and  77J  =  M  element.  Because  of  this  symmetry  and  the  relation 


e“  =  cos(:r)  +  i  sin(x)  (A. 7) 

the  odd  function,  isin(x),  will  have  no  contribution  when  the  displacement  is  calcu¬ 
lated.  The  final  simplified  matrix  are  presented  in  the  vector  form. 

The  middle  elements  are 

m=l:(Y  -  1),  n=l:(Y  -  1) 

C„,(n  -t-  Nm)  =  =  4  cos(K:„Xi)  cos(k„j/j)  (A.8) 

and  the  edge  elements  are 
m=0,  n=0 


m=0,  -  1) 


C,,(0)  =  =  1 


(A.9) 
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n  ^ 

m=0, 


^  n 
m=— ,  n=0 


C^(n)  =  =  26^"^^  =  2cos(K„yj) 


c„(f )  =  »  =  oos(«^w) 


C.(~)  =  =  e“¥*‘  =  cos(«^ii) 


m=— ,  n=l:(—  -  1) 


iVM. 


Cu,(n  +  -^-;r)  =  2e’'‘  *  =  2  cos(km Xj)  cos(«:„?/j) 

A  ^ 


M  N 
m=— ,  n=— 
2  ’  2 


.(y  +  jy)  =  =  cos(/CMXi)cos(K^7/j) 


,M 


m=l:(y  -  1),  n=0 


Cu,(— m)  =  26’''”*®*''"’®*'^  =  2e*''"*®’  =  2cos{KmXi) 
2 


,M  ,,  N 
m=l;(y  -  1),  n=y 


=  2e’''’"*‘''’*''‘^*'^  =  2 cos(KmXi)  cos(k ivy j) 

z  z 


(A.IO) 

(A.11) 

(A.12) 

(A.13) 

(A.14) 

(A.15) 

(A.16) 
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With  these  relations,  C,„  can  be  expressed  in  the  simplified  vector  form  in  Chap¬ 
ter  5.  For  m  =  0,  . . . ,  y  and  n  =  0,  . . . ,  f  ,  each  row  becomes 

Cu,  = 

[  1  •••  2cos(/c„j/j)  •••  cos{KN/2yj) 

2  cos(KmXi)  •  •  •  4cos(«;,„a:t)  cos{Knyj)  •  •  •  2  cos(/Cma:t)  cos(/CAr/2^j) 

COs{K,M/2Xi)  •••  2c0s{KM/2Xi)C0s{Knyj)  COs(KM/2Xi)  COs(KN/2yj)  ] 

(A.17) 


Appendix  B 

Existing  Formulation  for  Backward 
Filter 


This  computational  procedure  is  used  by  Bryson  [11]  to  solve  the  smoother  equa¬ 
tions. 

Continuing  from  the  measurement  downdate  relations  Eq.  6.16,  the  time  downdate 
procedure  starts  with  the  boundary  conditions 


S^{NT)  =  S{NT) 
X^{NT)  =  S(NT} 


(B.l) 


and  results  in  the  optimal  estimation  z(0)  at  n  =  0  using  the  boundary  condition 
Eq.  6.10  as 

z(0)  =  [S®"'  -b  S(0)]"'  [A^(0)  S(0)z(0)]  (B.2) 

With  the  boundary  conditions,  we  can  solve  the  two  point  boundary  value  problem 
by  sweeping  forward  from  n  =  0.  However,  the  Euler-Lagrange  equation  is  unstable 
in  either  direction.  A  stable  sequencing  can  be  done  by  storing  a  vector  f^{n)  and  a 
matrix  K®(n)  during  the  backward  sequencing 
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Kb(7i)  =  ■Wjl(NT-n)-Wj,(NT-n) 

(B.3) 

/^(n)  =  W7/(iVT-n)r^A^(n+l) 

The  forward  sequence  is  then 

/(n)  =  /®(ri)  -  K®(n)z(n) 

s(n)  =  Cz(n)  (B-4) 


z(n  +  1)  =  ^z(n)  +  r/(n) 


Appendix  C 

Rank- Two  Update  Procedure 


The  outer  loop  of  the  inverse  problem  solver  uses  the  Hessian  matrix  to  update 
the  estimate  of  the  impact  location.  The  impact  is  updated  using  the  quasi-newton 
method 

Xe(z  +  1)  =  Xe(i)  -  a3;^(i)Jx,(i)  (C.l) 

where  Xg  is  the  estimate  of  the  impact  location,  and  a  is  the  step  size  scaling  variable. 

The  Hessian  matrix,  jix(i),  is  calculated  approximately  using  a  rank-two  update 
procedure.  The  procedure  uses  the  gradiant  information  from  the  current  and  previous 
estimates  of  the  impact  location.  This  technique  is  shown  in  completion  in  Cuthbert 
and  Luenberger  [14]  and  [15]. 

The  change  in  the  estimated  impact  location,  p,  and  the  change  in  the  gradiant, 
g,  from  one  estimated  location  to  the  next  is 

p(i)  =  Xe(i  +  1)  -  Xe(i)  .^2) 

q(i)  =  jie(j  +  l)- jx.(*) 

The  update  uses  two  rank-one  matrices  which  are  the  outer  products  of  the  p  and  q 
matrices.  This  provides  at  most  a  rank-two  update. 


Jxx  (i  -t-  1)  ~  JizCO  "b 


(C.3) 


q^(i)p(i)  p(zrjg(i)p(i) 

The  initial  value  of  the  Hessian,  jii(O),  can  be  chosen  as  any  symmetric  posititve 
definite  matrix.  The  identitiy  matrix  is  commonly  used  causing  the  first  update  to 
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be  in  the  direction  of  steepest  descent. 

To  provide  for  easier  calculation  Eq.  C.3  can  be  inverted  directly.  The  approxi¬ 
mation  to  the  inverse  Hessian  matrix  H(z)  can  be  calculated  with 


H{i  -t- 1) 


^(^  4.  q^(0H(z)q(z)\  p(i)p^(i) 
MW +1^1+ 

pr(i)q(i)  (p(i)q"'(>)H(i)  +  H(i)q(.)p’'(i)) 


{C.4) 
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